home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 12 / Mac Magazin and MacEasy Magazine CD - Issue 12.iso / Sharewarebibliothek / Anwendungen / Wissenschaft & Technik / Yorick / yorick11-ppc folder / linpack.i next >
Text File  |  1995-04-17  |  7KB  |  255 lines
  1. /* Yorick equivalent of the Linpack benchmark */
  2.  
  3. func linpack {
  4.  
  5. elapsed= array(0.0, 3);
  6. old = second(0.0);
  7.  
  8. /* set the size of the matrix */
  9. if(is_void(is_a_mac)) {
  10.   n= 100;
  11. } else {
  12.   val= is_a_mac();
  13.   if(val == 1) {
  14.     /* run small problems on 68K Macs with no FPU */
  15.     n= 25;
  16.     write,"68K Macintosh without FPU";
  17.   } else if(val == 2) {
  18.     /* run mdeium problems on 68K Macs with FPU */
  19.     n= 50;
  20.     write,"68K Macintosh with FPU";
  21.   } else {
  22.     /* run full problem on Power Macs */
  23.     n= 100;
  24.     write,"Power Macintosh";
  25.   }
  26. }
  27. n= min(500, max(20, n));
  28. write," using size ",n;
  29. cray = .056;
  30. ops = (2.0*n^3)/3.0 + 2.0*n^2;
  31. /* WARNING WARNING!!! will not allocate extra storage,
  32.    AND will not test case where there is extra storage in matrix a */
  33. msiz= n;
  34.  
  35. /* allocate arrays */
  36. a= array(0.0, msiz, msiz);
  37. aa= array(0.0, msiz, msiz);
  38. b= array(0.0, msiz);
  39. bb= array(0.0, msiz);
  40. x= array(0.0, msiz);
  41. time= array(0.0, 8, 6);
  42. lda = msiz;
  43. norma= 0.0;
  44.  
  45. /* set the number of passes to make through the solver loop
  46.    so that time can be measured accurately
  47. */
  48. ntimes= 50;
  49.  
  50. /* generate matrix a and rhs b once and for all */
  51. matgen,aa,lda,n,bb,norma;
  52.  
  53. /* copy matrix into place and run first test */
  54. b= bb;
  55. a= aa;
  56. t1 = second(old);
  57. b= LUsolve(a,b);
  58. total = second(old) - t1;
  59.  
  60. /* compute a residual to verify results. */
  61.  
  62. x= b;
  63. a= aa;
  64. b= -bb;
  65. dmxpy,n,b,n,lda,x,a;
  66. resid= max(abs(b));
  67. normx= max(abs(x));
  68. eps = epslon(1.0);
  69. residn = resid/( n*norma*normx*eps );
  70. write,"    norm. resid      resid           machep",
  71.       "        x(1)          x(n)";
  72. write,format="%16.8e%16.8e%16.8e%16.8e%16.8e", residn,resid,eps,x(1),x(n);
  73.  
  74. write,format="\n\n    times are reported for matrices of order %6d\n",n;
  75. write,"     dgesv     mflops       unit      ratio";
  76.  
  77. time(1,3) = total;
  78. time(1,4) = ops/(1.0e6*total);
  79. time(1,5) = 2.0/time(1,4);
  80. time(1,6) = total/cray;
  81. write,format=" times for array with leading dimension of%6d\n", lda;
  82. write,format="%11.3e%11.3e%11.3e%11.3e\n", time(1,3),time(1,4),
  83.     time(1,5),time(1,6);
  84.  
  85. b= bb;
  86. a= aa;
  87. t1 = second(old);
  88. b= LUsolve(a,b);
  89. total = second(old) - t1;
  90. time(2,3) = total;
  91. time(2,4) = ops/(1.0e6*total);
  92. time(2,5) = 2.0/time(2,4);
  93. time(2,6) = total/cray;
  94.  
  95. b= bb;
  96. a= aa;
  97. t1 = second(old);
  98. b= LUsolve(a,b);
  99. total = second(old) - t1;
  100. time(3,3) = total;
  101. time(3,4) = ops/(1.0e6*total);
  102. time(3,5) = 2.0/time(3,4);
  103. time(3,6) = total/cray;
  104.  
  105. tm2 = 0;
  106. t1 = second(old);
  107. for(i = 1; i <= ntimes; i++) {
  108.   tm = second(old);
  109.   b= bb;
  110.   a= aa;
  111.   tm2 += second(old) - tm;
  112.   b= LUsolve(a,b);
  113. }
  114. total = (second(old) - t1 - tm2)/ntimes;
  115. total;total+tm2/ntimes;
  116. time(4,3) = total;
  117. time(4,4) = ops/(1.0e6*total);
  118. time(4,5) = 2.0/time(4,4);
  119. time(4,6) = total/cray;
  120.  
  121. write,format="%11.3e%11.3e%11.3e%11.3e\n", time(2,3),time(2,4),
  122.     time(2,5),time(2,6);
  123. write,format="%11.3e%11.3e%11.3e%11.3e\n", time(3,3),time(3,4),
  124.     time(3,5),time(3,6);
  125. write,format=" result of %4d passes through solver\n", ntimes;
  126. write,format="%11.3e%11.3e%11.3e%11.3e\n", time(4,3),time(4,4),
  127.     time(4,5),time(4,6);
  128. }
  129.  
  130. func second(old) {
  131.   extern elapsed;
  132.  
  133.   /* return the total CPU plus system time for this job,
  134.      less the input argument. */
  135.   timer,elapsed;
  136.   return elapsed(1)+elapsed(2)-old;
  137. }
  138.  
  139. func dmxpy (n1, y, n2, ldm, x, m)
  140. {
  141. /*   purpose:
  142.      multiply matrix m times vector x and add the result to vector y.
  143.  
  144.    parameters:
  145.  
  146.      n1 integer, number of elements in vector y, and number of rows in
  147.          matrix m
  148.  
  149.      y double precision(n1), vector of length n1 to which is added
  150.          the product m*x
  151.  
  152.      n2 integer, number of elements in vector x, and number of columns
  153.          in matrix m
  154.  
  155.      ldm integer, leading dimension of array m
  156.  
  157.      x double precision(n2), vector of length n2
  158.  
  159.      m double precision(ldm,n2), matrix of n1 rows and n2 columns
  160.  
  161.  ----------------------------------------------------------------------
  162. */
  163.   /*   cleanup odd vector */
  164.   j = n2%2;
  165.   if (j >= 1) {
  166.     y(1:n1) += x(j)*m(1:n1,j);
  167.   }
  168.   /*   cleanup odd group of two vectors */
  169.   j = n2%4;
  170.   if (j >= 2) {
  171.     y(1:n1) += x(j-1)*m(1:n1,j-1) + x(j)*m(1:n1,j);
  172.   }
  173.   /*   cleanup odd group of four vectors */
  174.   j = n2%8;
  175.   if (j >= 4) {
  176.     y(1:n1) += x(j-3)*m(1:n1,j-3)  + x(j-2)*m(1:n1,j-2) +
  177.                x(j-1)*m(1:n1,j-1) + x(j)*m(1:n1,j);
  178.   }
  179.   /*   cleanup odd group of eight vectors */
  180.   j = n2%16;
  181.   if (j >= 8) {
  182.     y(1:n1) += x(j-7)*m(1:n1,j-7) + x(j-6)*m(1:n1,j-6) +
  183.                x(j-5)*m(1:n1,j-5) + x(j-4)*m(1:n1,j-4) +
  184.                x(j-3)*m(1:n1,j-3) + x(j-2)*m(1:n1,j-2) +
  185.                x(j-1)*m(1:n1,j-1) + x(j)  *m(1:n1,j);
  186.   }
  187.   /*   main loop - groups of sixteen vectors */
  188.   jmin = j+16;
  189.   for(j = jmin; j <= n2; j += 16) {
  190.     y(1:n1) += x(j-15)*m(i,j-15) + x(j-14)*m(1:n1,j-14) +
  191.                x(j-13)*m(1:n1,j-13) + x(j-12)*m(1:n1,j-12) +
  192.                x(j-11)*m(1:n1,j-11) + x(j-10)*m(1:n1,j-10) +
  193.                x(j- 9)*m(1:n1,j- 9) + x(j- 8)*m(1:n1,j- 8) +
  194.                x(j- 7)*m(1:n1,j- 7) + x(j- 6)*m(1:n1,j- 6) +
  195.                x(j- 5)*m(1:n1,j- 5) + x(j- 4)*m(1:n1,j- 4) +
  196.                x(j- 3)*m(1:n1,j- 3) + x(j- 2)*m(1:n1,j- 2) +
  197.                x(j- 1)*m(1:n1,j- 1) + x(j)   *m(1:n1,j);
  198.   }
  199. }
  200.  
  201. func matgen(a,lda,n,b,&norma)
  202. {
  203.   init = 1325;
  204.   for(j = 1; j <= n; j++) {
  205.     for(i = 1; i <= n; i++) {
  206.       init = 3125*init % 65536;
  207.       a(i,j) = (init - 32768.0)/16384.0;
  208.     }
  209.   }
  210.   norma= max(a(1:n,1:n));
  211.   b(1:n) = a(1:n,sum:1:n);
  212. }
  213.  
  214. func epslon (x)
  215. {
  216. /*
  217.      estimate unit roundoff in quantities of size x.
  218.  
  219.       double precision a,b,c,eps
  220.  
  221.      this program should function properly on all systems
  222.      satisfying the following two assumptions,
  223.         1.  the base used in representing dfloating point
  224.             numbers is not a power of three.
  225.         2.  the quantity  a  in statement 10 is represented to
  226.             the accuracy used in dfloating point variables
  227.             that are stored in memory.
  228.      the statement number 10 and the go to 10 are intended to
  229.      force optimizing compilers to generate code satisfying
  230.      assumption 2.
  231.      under these assumptions, it should be true that,
  232.             a  is not exactly equal to four-thirds,
  233.             b  has a zero for its last bit or digit,
  234.             c  is not exactly equal to one,
  235.             eps  measures the separation of 1.0 from
  236.                  the next larger dfloating point number.
  237.      the developers of eispack would appreciate being informed
  238.      about any systems where these assumptions do not hold.
  239.  
  240.      *****************************************************************
  241.      this routine is one of the auxiliary routines used by eispack iii
  242.      to avoid machine dependencies.
  243.      *****************************************************************
  244.  
  245.      this version dated 4/6/83.
  246. */
  247.   a = 4.0/3.0;
  248.   do {
  249.     b = a - 1.0;
  250.     c = b + b + b;
  251.     eps = abs(c-1.0);
  252.   } while(eps == 0.0) ;
  253.   return eps*abs(x);
  254. }
  255.